NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


Technical Memorandum 33-724 

An Initial Feasibility Stage for Stoer’s 
Constrained Least Squares Algorithm 

Charles L. Lawson 


f {NASA-eB-142189) AN INITIAL FEASIBILITI 
STAGE FOR STOER'S CONSTRAINED LEAST SQUARES 
ALGORITHM (J[et Propalsion Lab.) 19 p HC . 

$3.25 CSCL 12A 

I G3/64 


N75-2202^ 


Qnclas 
19428 J 



JET PROPULSION LABORATORY 
CALtFORNIA INSTITUTE OF TECHNOLOGY 
PASADENA, CALIFORNIA 


April 15, 1975 




/ 


Prepared Under Contract No, NAS 7*100 
National Aeronautics and Space Administration 



PREFACE 


The work described in this report was performed by the Data Systems 
Division of the Jet Propulsion Laboratory. 


JPL Technical Memorandum 33-724 


ui 



ACKNOWLEDGMENTS 


The author has benefited from discussions with Prof, Bartels at the 
symposium cited in [l] and also from access to an Algol code for Stoer's 
algorithm written by Mr. P, Patzelt and kindly provided by Prof. Stoer. 


IV 


JPL Technical Memorandum 33-724 



CONTENTS 


1. Introduction 1 

2. Computing an initial feasible vector Xq 3 

3. Partial triangularization of N 6 

4. Initialization of Stoer’s state variables 11 

References * \z 


JPL Technical Memorandum 33-724 


V 



ABSTRACT 


A procedure is described for computing an initial feasible vector, Xq, 
for Steer's algorithm for solving the linear least squares problem subject to 
linear equality and inequality constraints. The procedure described fits well 
with Stoer's algorithm since much of the computation performed to determine 
Xq accomplishes initializing transformations of the problem data, which 
would otherwise be done in Stoer's algorithm after being given an Xq. 
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AN INITIAL FEASIBILITY STAGE FOR STOER'S CONSTRAINED 


LEAST SQUARES ALGORITHM 
CHARLES L. LAWSON 


I . Introduction. 

Denoting real matrices and vectors by upper and lower case letters 
respectively, introduce the following constrained least squares problems: 

(1) Problem LSIE Least squares with inequality and equality 
constraints, 

(la) Minimize |i Ex - f II 
subject to 

(lb) Cx = d. 

and 

( 1 c) Gx > h 

(2) Problem LSE Least squares with equality constraints 
[Eq (la) and (lb)j. 

(3) Problem LSI Least squares with inequality constraints 
[Eq (la) and (Ic)]. 

(4) Problem NNLS Nonnegative least squares. 

(4a) Minimize II Ex - fll 

subject to 

(4b) X > 0 
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(5) Problem L>DP Least distance programming. 


(6a) Minimize llxll 

subject to 

(5b) Gx > h 

Stoer [3], Lawson and Hanson [2, Chapter 23] and Bartels [l] have given 

algorithms for Problem LSIE which work directly with the matrix E rather 

T 

than with the nonnegative -definite symmetric matrix E E which is used in most 
quadratic programming algorithms. 

Although the algorithms described by Lawson and Hanson and by Bartels 
are different, they have the common feature that Problem LSIE is transformed 
via duality relations to a Problem NNLS. (Note that in this context the symbols 
E, X, and f of Eq (4a) do not denote the same quantities as in Eq (la) ). Prob- 
lem NNLS is then solved by a descent algorithm having finite convergence. 

One feature of these dual approaches is that there is no special problem 
of finding an initial feasible vector. The vector x - 0 is feasible for 
Problem NNLS. 

In contrast the Stoer algorithm maybe classified as a primal algorithm. 

It deals directly with the objective function of Problem LSIE, also by a 
descent algorithm having finite convergence. The Stoer algorithm requires an 
initial feasible vector, x^, satisfying the constraints of Eq (lb) and (Ic), 

Two possibilities are mentioned in [3] regarding the determination of x^. 
One is that a feasible vector may be known a priori. The other is that if the 
LSIE algorithm is extended to include a linear programming capability (such 

"‘'More precisely Bartels discusses the problem with generality which encom- 
passes either working with E or with E^^E, 
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an extension is included in loc. cit. ) then a special preliminary linear pro- 
gramming problem could be solved to determine Xq. This latter approach 
would apparently require the temporary introduction of some artificial vari- 
ables with an associated requirement for some additional two-dimensional 
arrays of working storage, however this storage question is not detailed in 
loc. cit. 

The purpose of this paper is to show how the LDP algorithm of [2] can be 
adapted for use as an initial feasibility stage for Stoer's algorithm. This 
approach requires only one additional row and/or column in some of the two- 
dimensional arrays already present in Stoer's algorithm. Furthermore it 
blends conveniently with the initialization phase of his algorithm in the follow- 
ing sense: The initialization phase tricingularizes the rows of the constraint 
equations (lb) and (Ic) which are satisfied as equations by the given x^. The 
proposed use of algorithm LDP will do most of the work of this triangularization 

as a by-product of determining x^. 

2. Computing an initial feasible v&ctor Xq. 

Assume the data denoted by E, f, C, d, G, and h in Eq (la)-(lc) is 
given, specifying a Problem LSIE. Assume the matrix C has m rows and 
n columns. To avoid complications in the exposition assume further that m < n 
and Rank (C) - m. 

Let K be an n X n orthogonal matrix that triangular! zes C from the right, 

thus 



where Lj^ is an m X m lower triangular nonsingular matrix. 
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With the change of variables 


( 7 ) X = Ky 

Eq (lb) and (Ic) are respectively transformed to 

( 8 ) [Lj : O] y = d 

and 

^ 9 ) [M : N] y 5 h 

Let y^ denote the first m components of y, and 72 denote the remaining n - m 
components. Equation ( 8 ) determines a unique value for yj, say y^, as the 
solution of 

( 10 ) Ljy^ = d 

Substituting y^ into Eq (9) gives 

(11) N y^ 2 : h - M y^ 

Define 

(12) p = h - M 

so that Eq (11 ) may be written as 

(13) Ny^>p 

We now seek an (n - m)-vector y^ satisfying Eq (13). We will impose 
the additional condition of seeking the y^ of least Euclidean norm satisfying 
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Eq (13). This is an LDP problem. As is shown in Chapter 23 of [2] the solu- 
tion to this problem can be extracted from the residual vector of a dual NNLS 
problem. 

The NNLS problem to be solved is 

(14) minimize 1|[0, 0, l] - v^[N : p]j| 

subject to 

(15) v > 0 

Letting v be the solution vector determined by Algorithm NNLS [2, p. 16l] 
for this problem, define the residual (n - m + l)-vector 

(16) = [0, 0, 1] - 0^[N : p] 

and let p denote the last component of Theoretically ^ will be zero if and 

only if II ^ It = 0 which will occur if and only if the constraint system of Eq (13) 
is inconsistent^ Otherwise with p ^ 0 define the (n - m)- vector y^ 

A y A 
= -r/p 

It can be shown [e. g, , Z, Chap. 23] that y is the minimal norm solution of 
Eq (13). Furthermore using the transformation of Eq (7) it follows that the 
n-ve ctor 

(18) ^ = K 



(17) 


A 

^2 

-1 


JPL Technical Memorandum 33-724 


5 



satisfies the constraints of Eq (lb) and (Ic) and thus is suitable for use as an 
initial feasible vector Xq in Stoer's algorithm. 

3. Partial triangularization of N. 

To initialize Stoer's algorithm we next wish to transform N to a partially 
triangularized form. Specifically let I denote the set of row indices of G for 
which X satisfies Eq (Ic) with equality. To initialize Stoer's algorithm, we 
must triangularize the rows of N with indices i€l. Note that I is also the set of 
row indices of N for which y^ satisfies Eq (13) with equality. 

The key observation is that generally most of the work needed to achieve 
this triangularization will already have been done by Algorithm NNLS . 

To provide the details supporting this observation we must summarize 
some properties of Algorithm NNLS. Algorithm NNLS induces a partition of 
the row indices of [N : p] into two sets, say i and its complement. Supposed 
contains ic indices and, by relabeling if necessary, assume these indices are 
1, 2, • • • , Partition [N : p] into its first k rows and the remaining rows as 

(19) [N : p] = 



Algorithm NNLS generates an (n-m + l)x(n-m + 1) orthogonal matrix U 
satisfying 


( 20 ) 


^1 

Pi 



0 

0 

^2 

Pz 

u = 


^2 

P2 

0 

1 


T 

1 

T 

®2 

O' 
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where is a k X k lower triangular nonsingular matrix. If the solution 
vector V produced by Algorithm NNLS is also partitioned into v^ and v^ where 
Vj denotes the first k components, then these subvectors satisfy 


(21) 


aT , T 

V, = Sj 

(22) 


aT - 

Vj > 0 

(23) 


aT „ 

= 0 

and satisfies 

> 


(24) 


Ni = P, 

and 



(25) 


72 = P2 


From these last two equations it follows that tcl. In other words all 
rows of [N : p] which have been triangular! zed by Algorithm NNLS are rows 
that needed to be triangularized to initialize Steer's algorithm. In addition 
however, y^ may "accidentally" satisfy some of the rows of Eq (25) with 
equality. If so these rows must eventually be brought into the triang;alariza- 
tion to complete the initialization. 

First however we attend to eliminating the effect of p^ in Eq (20) since 

p^ is extraneous for the triangularization of N^ that we need. The method will 

be similar to Row Removal Method 1 given in Chap. 27 of [2], Apply Givens 

rotations from the right to Eq (20) to transform the vectors s^ and s, to zero. 

This will necessarily transform c to 1 or -1 since the row vector [s^ , , (r] 

1 2 
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has -unit norm. These Givens rotations should operate on pairs of columns in 
the order (n - m + 1 , n - m), (n-m + 1 , n-m- 1 (n-m + 1 , 1 ), to 
preserve the zero structure of the first two blocks of the first block-row of the 
right-side matrix in Eq (20). 

Calling the product of these n - m rotations V this transforms Eq (20) 
to 



Pi 


^3 

0 

s' 

■"2 

P2 

UV = 

A 3 


‘2 

_ 0 

1 _ 


_0 

0 

± 1 _ 


where is a k X k lower triangular nonsingular matrix. Since 

[O • • • 0, l] UV = [O ♦ • ' 0, ±l] it follows that the last column of UV is 

[O ' • • 0, ±l] . Since UV is orthogonal this implies that the last row of UV 

must be [O • ’ • 0, ±l]. Thus UV is of the form 


(27) 



with W orthogonal. Substituting Eq (27) into Eq (26) we observe that 


(2 8 ) 



and 


(29) 




'^3 

0 ■ 


w = 



”2. 
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This gives our desired triangularization of N^. 
Define z by 

(30) = Wz 
Then from Eq (24) - (25) z satisfies 

(31) [L^ : O] z = 

and 

(32) [Ag : Bj] z 2 


We return now to consideration of the possibility (probably rare in 
practice) that I ^ If there are k indices in I and k > k this would be evi- 
denced by z satisfying k - k rows of Eq (32) with equality. By relabeling if 
necessary we assume these are the leading k - k rows of Eq (32). Let be 
an orthogonal matrix that brings these additional rows into the triangularization. 
Thus 


(33) 


■^3 

0 ' 

W = 

1 

0 ' 

.^3 

®3. 

/L 

I 

> 

\ 


where L^ is k X k lower triangular. In general could be singular however 
the nondegeneracy assumption which Stoer invokes [3, p„ 3 84] would imply 
that L^ is nonsingular. 

Define w by 

(34) z - w 
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Then w satisfies 


[L^ : O] w = 


[^4 : w > 


where p^ denotes the first k components of p and denotes the remaining 
subvector of p. 

Combining the various orthogonal transformation matrices we may define 


the n X n orthogonal matrix Q by 


U 0 ^ 


‘Q 0 


0 ± 


fK o~| r I 01 ri o") 

0 

ij Lo ij Lo uj 1_0 vj 


W 2 0 


0 0 1 


K 0*] 

0 

0 ij 


I 0 ^ IT 0 m 


W 0 0 W, 0 


0 0 ±10 0 1 ! 


This matrix Q triangularizes C and the first k rows of G as follows. 


0 0 


Q = Ml 0 




where 
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■M, 


M, 


is a partitioning of M of Eq (6) into its first k rows and the remaining rows. 


4. 


Initialization of Stoer's state variables. 


The Stoer algorithm is described in [3] in terms of transformations of a 
7-tuple denoted by |x, j; Q, L, G, R, h| where the dots are added here to pre- 
clude confusion with symbols used in the present paper. The first four of these 
items may be initialized as 

. A 

i ; = X 


J := I 


Q := Q 


and 


L : = 


Li 0 


Ml 


The items G, R, and h depend upon E and f of Eq (la) and would be computed 
as in [3] . 


JPL Technical Memorandum 33-724 


11 


References 


R, H, Bartels, Constrained Least Squares, Quadratic Programming, 
Complementary Pivot Programming and Duality, Proceedings of the 
Eighth Annual Symposium on the Interface of Computer Science and 
Statistics, Health Sciences Computing Facility, University of California, 
Los Angeles, Feb, 13-14, 1974. 

C, L. Lawson and R. J. Hanson, Solving Least Squares Problems, 
Prentice -Hall, 1974, 

J. Stoer, On the Numerical Solution of Constrained Least Squares 
Problems, SIAM J, Numer. Anal., Vol, 8, No. 2, 1971, pp. 382-411. 


JPL Technical Memorandum 33-724 

NASA - JPL - ComI, L.A., Calif. 



